Advection equation demo

using LinearAlgebra, Plots, Printf, SparseArrays

This implementation is again in Julia. Since tridiagonal matrices are so prevalent and important, Julia has a special data type for them. If you are using python or Matlab you want to use the spdiags command. For example, in python:

import numpy as np
from scipy.sparse import spdiags

m = 4
data = np.array([-np.ones(m), np.zeros(m) , np.ones(m)]);
diags = np.array([-1, 0, 1])
A = spdiags(data, diags, m, m)

A.toarray()  # just to see what it looks like

h = 0.01
m = convert(Int64,1/h)-1;
k = 0.1
T = 10.
S = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)) #increase dim by 1
a = 1.0;

S |> Array;

Periodic problem with forward Euler

For $a > 0$

$$ \begin{cases} u_t + a u_{x} = 0,\\ u(x,0) = \eta(x),\\ u(0,t) = u(1,t). \end{cases} $$

# Initial condition
η = x -> exp.(-20*(x .-1/2).^2)
# Initial condition chosen so that u(x,t) = sin(2*pi*(x - t)), if a = 1
# η = x -> sin.(2*pi*x)
# u = (x,t) -> sin.(2*pi*(x.-t))

#1 (generic function with 1 method)

Build MOL matrix

S = sparse(S) # Need to convert S to a new data type to allow new entries
S[1,end] = -1
S[end,1]  = 1
S |> Array

100×100 Matrix{Float64}:
  0.0   1.0   0.0   0.0   0.0   0.0  …   0.0   0.0   0.0   0.0   0.0  -1.0
 -1.0   0.0   1.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0  -1.0   0.0   1.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0  -1.0   0.0   1.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0  -1.0   0.0   1.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0  -1.0   0.0  …   0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0  -1.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0  …   0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  ⋮                             ⋮    ⋱         ⋮                      
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0  …   0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      1.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   1.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0  …  -1.0   0.0   1.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0  -1.0   0.0   1.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0  -1.0   0.0   1.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0  -1.0   0.0   1.0
  1.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0  -1.0   0.0

k = h^2
B = I - (a*k/(2h)*S)

100×100 SparseMatrixCSC{Float64, Int64} with 300 stored entries:
⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈
⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀
⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦

n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
# If solution is known
# plot!(x, u(x,t), xaxis = [0,1], yaxis = [-1,2],lw=1,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter)

fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
    t += k
    U = B*U
    if mod(i-1,tb) ≈ 0.0
        IJulia.clear_output(true)
        plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
        # If solution is known
        # plot!(x, u(x,t), xaxis = [0,1], yaxis = [-1,2],lw=1,label = @sprintf("u(x,t), t = %1.2f",t))
        plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
    end
end

Create a gif

plot()
anim = Animation()
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter)
frame(anim)

fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
    t += k
    U = B*U
    if mod(i-1,tb) ≈ 0.0
        plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
        plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter)
        frame(anim)
    end
end
gif(anim,"advection_periodic.gif")

┌ Info: Saved animation to 
│   fn = /Users/thomastrogdon/Dropbox (uwamath)/Teaching/586/amath-586-2022/notebooks/advection_periodic.gif
└ @ Plots /Users/thomastrogdon/.julia/packages/Plots/zozYv/src/animation.jl:114

Periodic problem with Lax-Friedrichs

$$ \begin{cases} u_t + a u_{x} = 0,\\ u(x,0) = \eta(x),\\ u(0,t) = u(1,t). \end{cases} $$

# Initial condition
η = x -> exp.(-20*(x .-1/2).^2)
# Initial condition chosen so that u(x,t) = sin(2*pi*(x - t)), if a = 1
#η = x -> sin.(2*pi*x)
#u = (x,t) -> sin.(2*pi*(x.-t))

#3 (generic function with 1 method)

h = 0.01
m = convert(Int64,1/h)-1;
T = 10.
S = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)) #increase dim by 1
a = 1.0;

S = sparse(S) # Need to convert A₀ to a new data type to allow new entries
S[1,end] = -1
S[end,1]  = 1
C₀ = sparse(Tridiagonal(fill(0.5,m),fill(0.0,m+1),fill(0.5,m)))
C₀[1,end] = .5
C₀[end,1] = .5
k = h # need this to avoid too much numerical dissipation
B = sparse(C₀ - (a*k/(2h))*S);

n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
# If solution is known
# plot!(x, u(x,t), xaxis = [0,1], yaxis = [-1,2],lw=1,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter)

fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
    t += k
    U = B*U
    if mod(i-1,tb) ≈ 0.0
        IJulia.clear_output(true)
        plot(x, U, xaxis = [0,1], yaxis = [0,1],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
        # If solution is known
        # plot!(x, u(x,t), xaxis = [0,1], yaxis = [-1,2],lw=1,label = @sprintf("u(x,t), t = %1.2f",t))
        plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
    end
end

Create a gif

plot()
anim = Animation()
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter)
frame(anim)

fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
    t += k
    U = B*U
    if mod(i-1,tb) ≈ 0.0
        plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
        plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter)
        frame(anim)
    end
end
gif(anim,"advection_periodic_LF.gif")

┌ Info: Saved animation to 
│   fn = /Users/thomastrogdon/Dropbox (uwamath)/Teaching/586/amath-586-2022/notebooks/advection_periodic_LF.gif
└ @ Plots /Users/thomastrogdon/.julia/packages/Plots/zozYv/src/animation.jl:114

Dirichlet BC with Lax-Friedrichs

$$ \begin{cases} u_t + a u_{x} = 0,\\ u(x,0) = \eta(x),\\ u(0,t) = g_0(t). \end{cases} $$

Consider the method of lines discretization:

$$ U'(t) = - \frac{a}{2h} A U(t) + \frac{a}{2h} g(t)$$

where

$$A = \begin{bmatrix} 0 & 1 \\ -1 & 0 & 1 \\ & -1 & 0 & 1 \\ & & \ddots & & \ddots\\ &&& -1 & 0 & 1 \\ &&& 1 & -4 & 3 \end{bmatrix}, \quad g(t) = \begin{bmatrix} g_0(t) \\ 0 \\ \vdots \\ 0 \end{bmatrix}.$$

Forward Euler produces:

$$U^{n+1} = U^n - \frac{ak}{2h} A U^n + \frac{ak}{2h} g(t). $$

We introduce a Lax-Friedrichs stabilization via a matrix $C$

$$ C =\begin{bmatrix} 0 & 1/2 \\ 1/2 & 0 & 1/2 \\ & 1/2 & 0 & 1/2 \\ & & \ddots & & \ddots\\ &&& 1/2 & 0 & 1/2 \\ &&& & 0 & 1 \end{bmatrix}.$$

Using the approximation $U^n \approx C U^n + g(t)/2$ we find the method

$$U^{n+1} = C U^n - \frac{ak}{2h} A U^n + \left( \frac 1 2 + \frac{ak}{2h}\right) g(t). $$

# Initial and boundary conditions
η = x -> exp.(-20*(x .-1/2).^2)
g0 = t -> sin(4*t)+1.

#7 (generic function with 1 method)

η(0.0) |> display
g0(0)

0.006737946999085467

1.0

h = 0.01
a = 1.0;
m = convert(Int64,1/h)-1;
k = h/(a)
T = 10.
A = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)); #increase dim by 1
A = sparse(A)
A[end,end-2:end] = [1,-4,3];

C = Tridiagonal(fill(0.5,m),fill(0.0,m+1),fill(0.5,m))
C = sparse(C)
C[end,:] *= 0.0
C[end,end] = 1.0
B = sparse(C - (a*k/(2h))*A);

n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[g0(t),U[end]], label = "BCs", seriestype = :scatter)

fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
    t += k
    U = B*U
    U[1] += (1/2 + a*k/(2h))*g0(t-k)
    if mod(i-1,tb) ≈ 0.0
        IJulia.clear_output(true)
        plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
        plot!([0,1],[g0(t),U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
    end
end

plot()
anim = Animation()
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[g0(t),U[end]], label = "BCs", seriestype = :scatter)
frame(anim)

fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
    t += k
    U = B*U
    U[1] += (1 + a*k/h)*g0(t)/2
    if mod(i-1,tb) ≈ 0.0
        plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
        plot!([0,1],[g0(t),U[end]], label = "BCs", seriestype = :scatter)
        frame(anim)
    end
end
gif(anim,"advection_LF.gif")

┌ Info: Saved animation to 
│   fn = /Users/thomastrogdon/Dropbox (uwamath)/Teaching/586/amath-586-2022/notebooks/advection_LF.gif
└ @ Plots /Users/thomastrogdon/.julia/packages/Plots/zozYv/src/animation.jl:114

Dirichlet BC with Lax-Wendroff

$$ \begin{cases} u_t + a u_{x} = 0,\\ u(x,0) = \eta(x),\\ u(0,t) = g_0(t). \end{cases} $$

# Initial and boundary conditions
η = x -> exp.(-40*(x .-1/2).^2)
g0 = t -> sin(4*t) + 1
# Initial condition chosen so that u(x,t) = sin(2*pi*(x - t)), if a = 1
# η = x -> sin.(2*pi*x)
# g0(t) = sin.(-2*pi*t)
# u = (x,t) -> sin.(2*pi*(x.-t))

#11 (generic function with 1 method)

h = 0.01
a = 1.0;
m = convert(Int64,1/h)-1;
k = h/(a)
T = 10.
A = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)) |> sparse; 
D = Tridiagonal(fill(1.0,m),fill(-2.0,m+1),fill(1.0,m)) |> sparse;
vec = [1,-4,3]
A[end,end-2:end] = vec; # same as Lax-Friedrichs
# Construct a backward second-order approximation of the second derivative
D[end,:] *= 0.0
D[end,end-3:end-1] += (vec[1]/2)*[-.5,0,.5]
D[end,end-2:end] += (vec[2]/2)*[-.5,0,.5]
D[end,end-2:end] += vec[3]*vec/4;

n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
# If solution is known
# plot!(x, u(x,t), xaxis = [0,1], yaxis = [-1,2],lw=1,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[g0(t),U[end]], label = "BCs", seriestype = :scatter)

fr = 1000 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
    t += k
    U = U - (a*k/(2h))*(A*U) + (a^2*k^2/(2h^2))*(D*U)
    U[1] += (a*k/(2h) + (a^2*k^2/(2h^2)))*g0(t-k)
    if mod(i-1,tb) ≈ 0.0
        IJulia.clear_output(true)
        plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
        # If solution is known
        # plot!(x, u(x,t), xaxis = [0,1], yaxis = [-1,2],lw=1,label = @sprintf("u(x,t), t = %1.2f",t))
        plot!([0,1],[g0(t),U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
    end
end

Convergence study with Lax-Wendroff

# Initial and boundary conditions chosen so that u(x,t) is the solution

F = x -> (exp.(-(x .- 1).^2) .+ 1).*sin.(2*pi*x)
a = 1.0
η = x -> F(x)
g0 = t -> F(-a*t)
u = (x,t) -> F(x.-a*t)

#19 (generic function with 1 method)

h = 0.2
p = 7

out = fill(0.0,p)
for i = 1:p
    h = h/2
    a = 1.0;
    m = convert(Int64,1/h)-1;
    k = h/(a)
    T = 10.
    A = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)) |> sparse; 
    D = Tridiagonal(fill(1.0,m),fill(-2.0,m+1),fill(1.0,m)) |> sparse;
    vec = [1,-4,3]
    A[end,end-2:end] = vec; # same as Lax-Friedrichs
    # Construct a backward second-order approximation of the second derivative
    D[end,:] *= 0.0
    D[end,end-3:end-1] += (vec[1]/2)*[-.5,0,.5]
    D[end,end-2:end] += (vec[2]/2)*[-.5,0,.5]
    D[end,end-2:end] += vec[3]*vec/4;
    
    n = convert(Int64,ceil(T/k))
    x = h:h:1 # include right end point
    U = η(x)
    t = 0.0

    for i = 2:n+1
        t += k
        U += -(a*k/(2h))*(A*U) + (a^2*k^2/(2h^2))*(D*U)
        U[1] += (a*k/(2h) + (a^2*k^2/(2h^2)))*g0(t-k)
    end
    out[i] = maximum(abs.(U - u(x,T)))
end
out[1:end-1]./out[2:end]

6-element Vector{Float64}:
 7.382313823819206
 7.841350172193508
 7.960061016385922
 7.989999195538282
 7.997501306853115
 7.999007598621097

η = x -> x.*exp.(-10*(x .- 1).^2)
η(0.00001)

η = x -> exp.(-20*(x .-1/2).^2)
g0 = t -> sin(4*t)^2.
h = 0.2

h = h/2^10
a = 1.0;
m = convert(Int64,1/h)-1;
k = h/(a)
T = 10.
A = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)) |> sparse; 
D = Tridiagonal(fill(1.0,m),fill(-2.0,m+1),fill(1.0,m)) |> sparse;
vec = [1,-4,3]
A[end,end-2:end] = vec; # same as Lax-Friedrichs
# Construct a backward second-order approximation of the second derivative
D[end,:] *= 0.0
D[end,end-3:end-1] += (vec[1]/2)*[-.5,0,.5]
D[end,end-2:end] += (vec[2]/2)*[-.5,0,.5]
D[end,end-2:end] += vec[3]*vec/4;
    
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0

for i = 2:n+1
    t += k
    U += -(a*k/(2h))*(A*U) + (a^2*k^2/(2h^2))*(D*U)
    U[1] += (a*k/(2h) + (a^2*k^2/(2h^2)))*g0(t-k)
end
U_true = U;

h = 0.2
p = 7
out = fill(0.0,p)
for i = 1:p
    h = h/2
    a = 1.0;
    m = convert(Int64,1/h)-1;
    k = h/(a)
    T = 10.
    A = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)) |> sparse; 
    D = Tridiagonal(fill(1.0,m),fill(-2.0,m+1),fill(1.0,m)) |> sparse;
    vec = [1,-4,3]
    A[end,end-2:end] = vec; # same as Lax-Friedrichs
    # Construct a backward second-order approximation of the second derivative
    D[end,:] *= 0.0
    D[end,end-3:end-1] += (vec[1]/2)*[-.5,0,.5]
    D[end,end-2:end] += (vec[2]/2)*[-.5,0,.5]
    D[end,end-2:end] += vec[3]*vec/4;
    
    n = convert(Int64,ceil(T/k))
    x = h:h:1 # include right end point
    U = η(x)
    t = 0.0

    for i = 2:n+1
        t += k
        U += -(a*k/(2h))*(A*U) + (a^2*k^2/(2h^2))*(D*U)
        U[1] += (a*k/(2h) + (a^2*k^2/(2h^2)))*g0(t-k)
    end
    U_true_p = U_true[end:-2^(10-i):1]
    out[i] = maximum(abs.(U[end:-1:1] - U_true_p))
end
out[1:end-1]./out[2:end]